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Abstract 

In a previous work we investigated the existence of Hopf degenerate 
bifurcation points for a differential delay equation modeling leukemia 
and we actually found Hopf points of codimension two for the consid- 
ered problem. If around such a point we vary two parameters (the 
considered problem has five parameters), then a Bautin bifurcation 
should occur. In this work we chose a Hopf point of codimension two 
for the considered problem and perform numerical integration for pa- 
rameters chosen in a neighborhood of the bifurcation point parameters. 
The results show that, indeed, we have a Bautin bifurcation in the cho- 
sen point. 
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1 Introduction 

The considered equation was taken from [9], [10] 
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and represents part of a model of periodic chronic myelogenous leukemia. 
The initial model consists of two delay equations, one for the density of 
proliferating cells, P, and one for the density of so-called "resting" cells, N. 
The latter equation is independent, hence it can be studied independently 
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(see [9], [IH])- Equation ([I]) is this second equation, with the unknown N 
made dimensionless by dividing it by a quantity with the same dimension. 
The parameters /3q, n, 5, k, r are positive real numbers. We do not insist 
here on their physical significance since this is largely presented in [9], |10| . 
Parameter k is of the form k = 2e _7r , with 7 positive. We take here, as in [6], 
k as an independent parameter, instead of 7. Note that, due to its definition, 
k < 2. We denote by a the vector of five parameters (a = (/3o, n, 5, k, r)). 
The equilibrium points of the problem are, as can be easily seen ([9], 

Xl = 0, x 2 = x 2 (a) = &k - 1) - l) 1 /™. 



The second one is acceptable from the biological point of view if and only if 

^(*-l)-l>0, (2) 

condition that implies k > 1. 

The equilibrium point X2 presents Hopf bifurcation for some points in 
the parameter space [9], [ID] . [5]. In [5] we developed the apparatus for 
investigating the normal form for a Hopf bifurcation point, by using the 
center manifold theory (we needed a second order approximation of the cen- 
ter manifold for this). In order to determine the normal form, we computed 
the first Lyapunov coefficient, l\(a). 

In [6] we searched for points of degenerate Hopf bifurcation, i.e. points 
a* in the parameter space with l\{a*) = 0. The method used relied also 
on the center manifold theory. We explored a quite extended zone of pa- 
rameters, having biological significance, and found that for n = 2 and each 
A) G {0.5, 1, 1.5, 2, 2.5}, k € {1.1, 1.2, 1.9}, values of r and 5 can be found, 
for which l\ = 0. We went further, and for these points we computed the 
second Lyapunov coefficient, by constructing a fourth order approximation 
of the center manifold. For all the points a* with l±(a*) = determined, 
we found < 0, and thus, by the definition in [TT], (x2, a*) represent 

Hopf points of codimension two. 

In the present paper we numerically investigate the occurrence of Bautin 
bifurcation in one of the Hopf points of codimension two found. In [6J 
we considered the restriction of the problem to a two-dimensional center 
manifold, this restriction is a two-dimensional problem and the theory of [7] 
works for its study. Some ideas concerning the restriction of the problem to 
the center manifold are presented here, in Section 2. 

In Section 3 we describe the Bautin bifurcation for two dimensional dy- 
namical systems, as it is presented in [7]. 

Then we chose a point in the parameter space for which l\ = and we 
explore its neighborhood in order to see whether we regain the bifurcation 
diagram of the Bautin bifurcation for our problem. For this, we numerically 
integrate our problem for the chosen parameters. As the graphs of the 
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trajectories obtained by numerical integration show, we have indeed a Bautin 
bifurcation in the chosen Hopf point of codimension two (Section 4). 



2 The problem restricted to the center manifold 

In [6] we considered the nontrivial equilibrium point %2 and the linearized 
equation around this point, that is (see also [3], [10J ) 

z(t) = —[Bi + S]z(t) + kBiz(t — r), (3) 

where z = x — X2{a), B\ = /3'(x2)x2 + /3(x2), and (3(x) = — — — . 

1 + x n 

The characteristic equation corresponding to (|3j) is 

A + 5 + B x = kB ie - Xr . (4) 

The eigenvalues depend on the vector of parameters, a. 

Assume that we have a point a* in the parameters space such that, for 
a in a neighborhood U of a*, there are two eigenvalues, that we denote by 
A a i,2 with the property that all other eigenvalues have negative real part, 
and, at a = a*, X a *,i,2 = 

The analysis in [5j shows that a* = (/3q, n*, 5* , k* , r*) satisfies the above 
condition if and only if the relation 

= a vc C os((5*+B*)/(k*B*)) 

y/(k*Bi) 2 - (5* +B x y ' 1 J 

is satisfied, where B* is the value of B\ at a*. 

For a = a* , a two-dimensional local invariant manifold (the local center 
manifold) exists and the reduction of the problem to this manifold leads to 
the ordinary differential equation: 

^=u*iu+ J2 i^9ik(a*)v?rf, (6) 

j+k>2 J ' 

where u : R i— > C. 

The formalism for the construction of an approximation of the center 
manifold and that for computing the coefficients gjk{o*) are fully presented 
in [6]. We remind here only some elements of that construction, that relies 
on the general ideas in [I], [2]. 

We considered the Banach space 

B = {ip : [— r, 0] i->- R, ij) is continuous on [— r, 0]} , 

and its complexification, denoted Be- We denoted by M. the subspace of 
Be spanned by the two eigenfunctions ip a * 1,2(3) = e ±UJ * ts , s G [— r, 0], corre- 
sponding to the two eigenvalues X a * 1,2 and by V a projector defined on Be, 
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with values in M. The local center manifold is locally invariant, tangent to 
M. in 0. It is the graph of a smooth function w a * : U C M. i-> (I-V)Bc, (U 
is a neighborhood of in M) that satisfies w a *(0) = 0. Thus a "point" <fi on 
the center manifold has the form cj) = zip a * i + z(p a * 2 + w a * (z(p a * 1 + ~z(p a * 2). 

For an initial condition on the center manifold, the solution x(-) of 
equation (P) satisfies 

X t = u(t)(p a * 1 + u(t)(f a * 2 + Wa* (u(t)<Pa* 1 + u{t)(p a * 2 ), 

where cct € is defined by xt(s) = x(t + s), s G [— r, 0], u(-) is the solution 
of equation © with the initial condition n(0) = no, and = uo^i + ^0922- 
When a £ U, the two eigenvalues A^i^ = ± iw(a) may have 

positive or negative real part. In each of these situations, there still is 
a two-dimensional local invariant manifold, a local unstable manifold when 
ReX a i,2 > 0, and a submanifold of the local stable manifold for Re\ a 1^ < 0. 
This latter case can be argued with the ideas of [8], adapted to the more 
simple case considered by us. Hence the solution of our problem has in this 
case also a representation of the form 

X t = U(t)(fal +U(t)<p a 2 + W a (u(t)(p a i +u(t)<p a2 ), 

where u satisfies an equation of the form 

= + ui)u + ^2 -^gajk(a)u J u k , (7) 

and w a is the function whose graph is the local invariant manifold. 

The above considerations show that the problem © (respectively ([7])) 
presents, for some initial value, a periodic solution iff the corresponding 
solution x(t) of (HJ) is periodic. Also, a solution of © (respectively ([7|)) 
spirals towards (or from 0), iff the corresponding solution of (JTJ) spirals 
towards (respectively from 0). Hence the study of equations ([6]), from the 
point of view of the Hopf or Bautin bifurcation leads to complete conclusions 
concerning these bifurcations for the problem ([1]). 

3 Bautin bifurcation for planar systems [7] 

Consider a system of two ODEs, that can be written as a single complex 
equation as 

z = X(a)z+ -^9jk(a)z j z k , (8) 

j+k>2^' 

where a = (01,02) £ M 2 . 

In the hypotheses that a certain value oq of a exists such that: 

• A(a ) = i^o, 
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• h(ao) = 0, 

• h(a ) / 0, 

• the map a — > (pi(a), p 2 (a)), where p\{a) = ^fey, P2{ot) = h(&) is 
regular at ao, 

equation ([8]) may brought by several transform of functions and of parame- 
ters to the form: 

u = (/Ui(a) + i)u + p 2 (a)u\u\ 2 + L 2 (a)u\u\ 4 + 0(\u\ 5 ), (9) 

where u : M. h-> C and L 2 (p) = l2(a(fj,)). 

Moreover, in [7J is proved that eq. ([9]) is locally topologically equivalent 
with 

ii = (b\ + i)u + b 2 u\u\ 2 + su|n| 4 , (10) 

where b\ = yUi, 62 = ^1-^2(^)1^2 and s is the signature of ^2(00)- 

In order to describe the phase portrait for the parameters varying around 
the point «o (equivalent to (b±, 62) = (0, 0)) it is useful to consider the polar 
form of the above equation: 

p = p(b\ + b 2 p 2 + sp 4 ), 
9 = 1. 

The limit cycles are obtained by solving the equation: 

h + b 2 p 2 + p 4 = 0, 

and, by studying the number of its solutions as function of b\,b 2 , the bifur- 
cation diagram of the Bautin bifurcation is obtained. 

We reproduce in Fig. 1 the bifurcation diagram for the case s = — 1 since 
for all our Bautin bifurcation points found the second Lyapunov coefficient 
is negative. 

We see that in a neighborhood of the origin in the plane of the parameters 
(61, 62) the phase portrait in a neighborhood of u = has very different 
aspects. These are described in [7], but for the sake of completeness, we 
point out a few ideas here. In the zone 1 of the Bautin bifurcation diagram, 
that lies between the b 2 < part of the axis 61 = and the curve T, the 
point u = is an attractive focus; when we cross the axis b\ = entering 
in the zone 2 (the half-plane b\ > 0) , a supercritical Hopf bifurcation takes 
place and a stable limit cycle occurs, while the point u = loses stability; 
then, starting from the first quadrant, when we cross the axis b\ = to 
arrive in the zone 3 (lying between the b 2 > part of the axis 61 = and 
the curve T), a subcritical Hopf bifurcation takes place, and an unstable 
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Figure 1: Bautin bifurcation diagram for sign^(ao) < 0. 



limit cycle is born, in the interior of that previously formed. Then, for 
the parameters (b\, 62) on the curve T, the two cycle collide in a single 
cycle, that is repulsive on its interior side and attractive on its exterior side, 
and after crossing the curve T, arriving in the zone 1 again, the two cycles 
disappear. 

Hence, the most interesting feature of this bifurcation is the presence of 
two limit cycles one inside the other for the parameters lying between the 
axis bi = and the curve T. The exterior cycle is stable (attractive) while 
the interior one is unstable (repulsive). 

3.1 Numerical confirmation of Bautin bifurcation 

We have chosen the case n* = 2, /3q = 2.5, k* = 1.01, that is not among 
those found in [6] . We took this case hoping to have a small r at the Bautin 
bifurcation and intending to have the second Lyapunov coefficient not very 
close to zero (this choice is justified by the table contained in Fig. 6 of [B]). 

By using the methods presented in [6], we find that at 
r* = 5.301432998, 6% = 0.0023073665, we have h = 0, while l 2 = -0.0662. 

In order to see if we can regain the bifurcation diagram above for our 
problem, we performed numerical integrations of the delay differential equa- 
tion (pQ) for the above values of n, (3q, k and values of r, 5 around r*, 5q. We 
used the routine dde23 of Matlab. 

In Fig. 2 we see, in the plane (r, S), the curve of points where u = 0, 
the point B where l\ = 0, and the points where we performed the numerical 
integration. 

In order to put into light the behavior of the solution, that is qualitatively 
described in the bifurcation diagram, for a chosen point in the parameters 
space, we have to take several initial functions situated at different distances 
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Figure 2: Curve of points (S, r) where ReX\^ = 0, for n = 2, /3 = 2.5, k = 1.01, 
the point of Bautin bifurcation (point B), and the points chosen for numerical 
integration. 




Figure 3: Phase portrait for the parameters in the point Pi, for c = 0.5. 



from the equilibrium point. We took the initial function (p of the form 
(p = X2 + ce^ s cos(us), where A a x,2 = A 4 ± iu; are the eigenvalues of the 
linearized problem at the chosen parameters, and c is a new parameter, that 
we vary. 

The results of the integrations, for each of the considered points and for 
some choices of c are represented vs time, but also in x(t), x(t) plots. 

We consider first the solutions for two points P±, P{ in the zone of the 
(<5, r) plane, corresponding to the zone 1 of the bifurcation diagram. More 
precisely the point Pi has the coordinates 5 = 0.002, r = 5.93, while for 
P{, 8 = 0.0024, r = 5.14 (see Fig. 2). 

The behavior of the solution for these two points is that corresponding 
to an attracting focus. Since we obtained qualitatively the same image for 
several values of c, we represent only one of these, that for c = 0.5, in Figs. 
3 and 4. 
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Figure 5: Phase portrait for the parameters in the point P2, for c = 0.001. 

The two points considered next, P2 (5 = 0.0024, r = 5.2) and P' 2 (5 = 
0.0015, r = 7.56) are situated in the zone corresponding to zone 2 of the 
bifurcation diagram (see Fig. 2). 

We see in Figs. 5-6 that a stable limit cycle occurs by Hopf bifurcation. 
In these two figures we present the behavior of the solution corresponding 
to P2 and c = 0.001, respectively c = 0.2. In Figs. 7 and 8 we present 
the behavior of the solution corresponding to P 2 and c = 0.1, respectively 
c = 5. We remark that at each oscillation, on this limit cycle, at the end of 
a descending branch, in the x versus t representations, a small superposed 
oscillation occurs, that produces a little spiral in the left of the x versus x 
representation. The moment when the solutions enters on this limit cycle 
depends on the distance between the initial function and xi- 

It was very difficult to find a point in the parameter plane, presenting 
the behavior corresponding to that of zone 3 of the bifurcation diagram. 
We suppose this is so because of the curve T that is deformed and may 
be very close to the curve oj = 0. However, we found that for the point 
P3 (d = 0.0015, r = 7.55), the behavior of the solution is the following: for 
initial functions close to X2, that is c < 0.42, the solution spirals towards 
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Figure 6: Phase portrait for the parameters in the point P2, for c = 0.2 
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Figure 7: Phase portrait for the parameters in the point P' 2l for c = 0.1 
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Figure 9: Phase portrait for the parameters in the point P3, for c = 0.1. 
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Figure 10: Phase portrait for the parameters in the point P3, for c = 0.42. 

X2, while for c > 0.425 it is visible that the solution spirals away from 
a repulsive limit cycle, having increasing (with time) amplitude. When 
time increases, the solution tends to a large attractive limit cycle (that 
was previously formed by Hopf bifurcation when passing from zone 1 to 
zone 2). This types of behavior are shown in Figs. 9-13 where we took 
c = o.l, c = 0.42, c = 0.425, c = 0.45, c = 0.6. In Figs. 10 and 11 the 
repulsive limit cycle is clearly visible, while in Figs. 12 and 13 the exterior 
attractive cycle is present. 

Hence, by numerical integrations we confirmed that the Bautin bifurca- 
tion, predicted by theoretical considerations, actually takes place, for the 
considered differential delay equation, in one of the points with l\ = 0, I2 < 
0. 

It is important to remark that a consequence of the Bautin bifurcation 
is, for a certain zone in the parameter space, the occurrence of a repulsive 
limit cycle inside an attractive one. There, the behavior of the solution of 
eq. (1) strongly depends on the initial function. If the initial function is 
close to the equilibrium point X2, the solution spirals towards x-i , while if 
the initial function is "far" from X2, it will spiral towards the exterior stable 
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Figure 11: Phase portrait for the parameters in the point P3, for c = 0.425. 
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Figure 12: Phase portrait for the parameters in the point P3, for c — 0.45. 
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Figure 13: Phase portrait for the parameters in the point P3, for c = 0.6. 
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limit cycle. From the point of view of the studied illness these two behaviors 
are quite different and this shows the importance of being able to control 
the initial condition. 
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